Response of electrically coupled spiking neurons: a cellular automaton approach 
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Experimental data suggest that some classes of spiking neurons in the first layers of sensory 
systems are electrically coupled via gap junctions or ephaptic interactions. When the electrical 
coupling is removed, the response function (firing rate vs. stimulus intensity) of the uncoupled 
neurons typically shows a decrease in dynamic range and sensitivity. In order to assess the effect 
of electrical coupling in the sensory periphery, we calculate the response to a Poisson stimulus 
of a chain of excitable neurons modeled by n-state Greenberg-Hastings cellular automata in two 
approximation levels. The single-site mean field approximation is shown to give poor results, failing 
to predict the absorbing state of the lattice, while the results for the pair approximation are in good 
agreement with computer simulations in the whole stimulus range. In particular, the dynamic range 
is substantially enlarged due to the propagation of excitable waves, which suggests a functional role 
for lateral electrical coupling. For probabilistic spike propagation the Hill exponent of the response 
function is a = 1, while for deterministic spike propagation we obtain a = 1/2, which is close to 
the experimental values of the psychophysical Stevens exponents for odor and light intensities. Our 
calculations are in qualitative agreement with experimental response functions of ganglion cells in 
the mammalian retina. 
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I. INTRODUCTION 

Unveiling how neuronal activity represents and pro- 
cesses sensory information remains a very difficult prob- 
lem, despite theoretical and experimental efforts under- 
taken by neuroscientists for the last several decades (for 
a recent review, see pj). In this broad context, rela- 
tively little attention has been devoted to the question 
of how organisms cope with the several orders of mag- 
nitude spanned by the intensities of sensory stimuli [2|. 
This astonishing ability is most easily revealed in humans 
by classical results in psychophysics 0] : the perception of 
the intensity of a given stimulus is experimentally shown 
to depend on the stimulus intensity ras~ log(r) (Weber- 
Fechner law) or ~ r a (Stevens law), where the Stevens 
exponent a is typically < 1. Those laws have in com- 
mon the fact that they are response functions with broad 
dynamic range, i.e. they map several decades of stimuli 
into a single decade of response. 

One would like to understand how this broad dynamic 
range is physically achieved by neuron assemblies. Re- 
cent experimental evidence suggests that electrical cou- 
pling among neurons in the early layers of sensory sys- 
tems plays an essential role in weak stimulus detec- 
tion. Deans et al. 4| showed that electrical coupling is 
present in the mammalian retina via gap junctions (ionic 
channels that connect neighboring cells). In particular, 
the spiking response of ganglion cells to light stimulus 
changes dramatically when the gap junctions are genet- 
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ically knocked out: both sensitivity and dynamic range 
are reduced 0. 

Another example comes from the olfactory system. 
The spiking response of isolated olfactory sensory neu- 
rons (OSNs) to varying odorant concentration usually 
presents a narrow dynamic range 0, Q . This is in con- 
trast with the response observed in the next layers of 
the olfactory bulb: both the glomerular 0,0 and mitral 
cell responses present a broader dynamic range than 
the OSNs. In this case, the tightly packed unmyelinated 
axons of OSNs in the olfactory nerve are believed to in- 
teract electrically via ephaptic interactions 01 (i- e - me- 
diated by current flow through the extracellular space), 
as shown by Bokil et al. [Tl|. In particular, their results 
indicate that a spike in a single axon can evoke spikes 
in all other axons of the bundle, suggesting that some 
computation is performed prior to the glomerular layer. 

Motivated by these results, previous papers have 
shown through numerical simulations that electrical cou- 
pling among neurons indeed changes the response func- 
tion in a way that is consistent with experimental re- 
sults. Due to the coupling, stimuli generate excitable 
waves which propagate through the neuron population. 
The interplay between wave creation and wave annihi- 
lation leads to a nonlinear amplification of the spiking 
response, increasing the sensitivity at low input levels 
and enhancing the dynamic range 0, In one di- 

mension, the robustness of the mechanism is attested 
by the diversity of models employed: cither the bio- 
physically realistic Hodgkin-Huxley equations [T^ . Il4| . 
a lattice of nonlinear coupled maps 0, 0, ^| , or the 
Greenberg-Hastings cellular automata (GHCA) [l2L IT^ 
yield qualitatively similar results. The same phenomenon 
has recently been observed in simulations with the two- 
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dimensional GHCA [P^ . 

In this paper we calculate the response of excitable 
GHCA model neurons 0] , where the bidirectional (elec- 
trical) coupling is modeled by a probability p of spike 
transmission. While the uncoupled case p — can be ex- 
actly solved, the coupled case p > is handled within two 
mean field approximations, namely at the single-site and 
pair levels. The aim is to shed light on the analytical be- 
havior of the response function for the one-dimensional 
case, therefore building on previous efforts which have 
relied entirely on numerical simulations. 

Our focus on the response of a continuously driven spa- 
tially extended excitable system should be carefully con- 
fronted with other recent studies, where the main inter- 
est has been on phase transitions between an excitable 
and a self-sustained collective state. For instance, the 
SIRS model of epidemics in hypercubic lattices has been 
recently investigated under the mean field and pair ap- 
proximations In those contagion models, station- 
ary self-sustained activity becomes stable for sufficiently 
strong connection among neighbors, a behavior which has 
been shown to be universal under very general assump- 
tions |20j| . Similar results have been obtained for a vari- 
ety of neuronal models, including collective responses to 
a localized transient stimulus |2lLl22| . as well as the emer- 
gence of sustained activity in complex networks [2^, |23| . 

While interesting in its own, the framework of stable- 
unstable collective transitions does not seem particularly 
suited for our modeling purposes. To account for sen- 
sory responses, the employed GHCA model is an ex- 
citable system which always returns to its absorbing state 
in the absence of stimulus, there are no phase transi- 
tions. The refractory period of the GHCA model neu- 
rons is absolute (unlike, say, reaction-diffusion lattices), 
mimicking the deterministic behavior of continuous-time 
systems like the Hodgkin-Huxley equations or integrate- 
and-fire models ^4[. The only source of stochasticity 
of the model regards the firing of the neurons. Stimuli 
can come from spiking neighbors (with probability p) or 
from an "external" source, which is modeled by a Poisson 
process and represents sensory input. Therefore, in the 
limit p = 1 the dynamics is that of a deterministic ex- 
citable lattice being stochastically stimulated, which casts 
the problem into the framework of probabilistic cellular 
automata [24| . 

The paper is organized as follows. In section [HI the 
GHCA rules are described; section IIIII contains the ex- 
act calculations for the response of uncoupled neurons, 
while in sections [Tvl and IVl results for the coupled case 
are discussed in the mean field and pair approximations, 
respectively. Our concluding remarks are presented in 
section IVI1 



II. THE MODEL 

In the n-state GHCA model for excitable sys- 
tems, the instantaneous membrane potential of the i-th 



cell (i = 1,...,L) at discrete time t is represented by 
Xi(t) € {0,1,..., n - 1}, n > 3. The state Xi(t) = 
denotes a neuron at its resting (polarized) potential, 
Xi(t) = 1 represents a spiking (depolarizing) neuron and 
Xi(t) — 2, ... ,n — I account for the afterspike refrac- 
tory period (hyperpolarization) . We employ the simplest 
rules of the automaton: if Xi(t) — 0, then Xi(t + 1) = 1 
only if there is a supra-threshold stimulus at site i; oth- 
erwise, Xi(t + 1) = 0. If Xi(t) > 1, then Xi(t + 1) = 
(xi(t) + l)mod n, regardless of the stimulus. In other 
words, the rules state that a neuron only spikes if stim- 
ulated, after which it undergoes an absolute refractory 
period before returning to rest. 

Whether the neurons are isolated or coupled is implicit 
in the definition of the supra-threshold stimulus. We as- 
sume external supra-threshold stimuli to be a Poisson 
process with rate r (events per second). Hence at each 
time step an external stimulus arrives with probability 

A(r) = 1 - (1) 

per neuron. Notice that t — 1 ms corresponds to the 
approximate duration of a spike and is the time scale 
adopted for the time step of the model. The number of 
states n therefore controls the duration of the refractory 
period (which corresponds to n — 2, in ms). In the bio- 
logical context, r could be related for example with the 
concentration of a given odorant presented to an olfac- 
tory epithelium 5], or the light intensity stimulating a 
retina |^y. We shall refer to r as the stimulus rate or 
intensity. 

When electrically coupled, neurons at rest can also be 
stimulated by their neighbors. We define p and q as 
the probabilities that a resting neuron spikes as a conse- 
quence of transmission (ionic current flow) from respec- 
tively one or two spiking neighbors [see Eq. . We keep 
p and q as two independent parameters in most calcula- 
tions to show the robustness of some asymptotic results. 
In the simulations, we concentrate on the more physically 
intuitive choice of g = 1 — (1— p) 2 , where the contribu- 
tions from two spiking neighbors are independent. 

Let pf^ (k) be the probability that the i-th neuron is 
in state k at time t. Since the dynamics of the refractory 
state is deterministic, the equations for k > 2 are simply 

P t « (2) = Pf (1) 
P t «(3) = P t W (2) 

P^n-l) = P t [l \n-2). (2) 

To describe the coupling among first neighbors, let 
P t (k, l,m) be the joint probability that sites i — 1, i 
and i + 1 are respectively in the states k, I and m at 
time t. Following the definitions of A, p and q above, the 
equation for Pj (1) thus becomes 
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where X x satisfies 



P t %(l) = [l-(l-A)(l-g)]P t w (1,0,1) 



(i), 



'n-1 



-[l-(l-A)(l-p)] £P t (i) (l,0,fc) 



E^ w ( fc >°,i) 



n— 1 n— 1 



+aEE p *°(w) 



(3) 



Finally, the dynamics for (0) can be obtained by the 
normalization condition 



Ei 5 t (i) (fe) = i,vt,i, 



(4) 



fe=0 



which completes the set of equations for one-site proba- 
bilities. 

It is reasonable to assume homogeneity in the sys- 
tem when L — > oo, so we can drop the superscript 
(i) in Eqs. (121 H and in what follows. We also ex- 
pect isotropy (right-left symmetry) in the probabilities: 
P t (k,l) =P t (l,k), Pt(k,l,m) = P t (m,l,k) etc. Recalling 
the normalization condition 2ji=o Pt{ji, hi • • ■ > 3 m) = 
Pt{j2, ■ ■ ■ ,jm), one can rewrite Eq. © as 



P t+1 (l) = APt(0) + 2p(l-A)P t (l,0) 

+ (l-A)(g-2p)P 4 (l,0,l). (5) 

The stationary value of any joint probability will 
be denoted by omitting the subscript i, thus P(») = 
limt^oo P* (•)■ We start by solving Eqs. ancl m 
the stationary state, which together yield 



P(A a 



xFrr 



(8) 



The dynamic range is therefore the number of decibels 
of input which are mapped into the ~ 9.5 dB of output 
comprised in the [0.1P max , 0.9F max ] interval (see Fig- EJt- 
In the biological context of the model, it measures the 
ability of the system to discriminate different orders of 
magnitude of stimulus intensity. We will show below that 
if one chooses to calculate 6 r using r x = — t _1 ln(l — A x ) 
instead of A x in Eq. results are essentially unchanged. 



III. UNCOUPLED NEURONS 

The uncoupled case p — q = can be exactly solved by 
taking the limit t — > oo in Eq. 101 which, together with 
Eq. ©, yields 



P(1) = /(A) = 



A 



l + (n-l)A 



(9) 



This linear saturating response is depicted for n = 3 
(Figs. EEl and n = 10 (Figs. UEJ), in complete agree- 
ment with simulations. It belongs to the family of Hill 
functions defined by H a (x) = Cx a / (xq + x a ), where the 
Hill exponent in this case is a = 1. 

The dynamic range can be promptly calculated: 
6\(n) = 101og 10 {[l + 9n]/[l + n/9]} and S r (n) = 
101og 10 {ln[l + 9/n]/ln[l + l/(9n)]}, both of which 
rapidly converge to 101og 10 (81) ~ 19 dB for moderate 
values of n (see lower curves in Fig.0J). As we shall see, 
the electrical coupling can lead to dynamic ranges typi- 
cally twice as large. 



P(0) = l-(n-l)P(l), (6) 

a result which is exact and holds Vp, q. 

We are interested in obtaining the behavior of P(l) as 
a function of A (or r). Note that P(l) coincides with the 
average firing rate per neuron (measured in spikes per ms, 
according to the choice of r) in the limit L,t — * oo. In 
simulations, firing rates have been calculated by division 
of the total number of spikes in the chain by LT, where 
T ~ O(10 5 ) and L ~ O(10 5 ) were the typical number of 
time steps and model neurons employed [l2j. We define 
P(A) = P(l) as the response function of the system. 

Due to the absolute nature of the refractory period, the 
maximum firing rate of the model neurons is F max = 1/n, 
a result which is easily obtained Vp, q by setting A = 1 in 
Eqs. JSJ and ©. The dynamic range 5\ of the response 
curve P(A) follows the definition commonly employed in 
biology [1 HI: 

*A = 101og 10 g£) , (7) 



IV. COUPLED NEURONS: MEAN FIELD 
APPROXIMATION 

As can be seen in Eq. (JSJ, Pt(l) depends on two- and 
three-site probabilities, and in general fc-site probabilities 
depend on up to (fc + 2)-site probabilities. The dynamical 
description of the system thus requires an infinite hierar- 
chy of equations. The mean field approximation at the 
single-site level corresponds to the simplest truncation of 
this hierarchy, and consists in discarding the influence of 
all neighbors in the conditional probabilities 26], thus 
Pt(h\j2, ■ ■ ■ ,3m) ~ Pt(ji), which leads to 

771 

Pt(jl,...,3m)~ Y[Pt(3k) ■ (10) 

fc=i 

In this approximation, Eq. JSJ) becomes 

P i+1 (l) w P t (0){A + 2p(l-A)P i (l) 

+ ( 9 -2p)(l-A)P t (l) 2 } , (11) 
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FIG. 1: Response curves for (a) n = 3 and (b) n — 10 au- 
tomata: simulations (symbols) and mean field approximation 
[lines, according to Eq. lilt ]. From bottom to top, p — 0, 
0.3, 0.6 and 1, q = 1 — (1 — p)' 2 . In the simulations, standard 
deviations over 10 runs are smaller than symbol sizes, so error 
bars are omitted in all figures. Notice the negative slope and 
multi-valuedness of the single-site approximation for p > 1/2 
and A < 0. 



which, together with Eq. JSJ, can be used to eliminate 
P(0) and render P(l) = F(X) implicitly through the re- 
lation 

(1 - 2p)F + (2pn - q)F 2 + (n - l)(g- 2p)F 3 
~ [1 - (n - \)F] [1 - 2pF + (2p - q)F 2 } ' 

(12) 

As a consistency check, notice that setting p = q = in 
Eq. (|12JI recovers Eq. (JHJ (in other words, mean field is 
exact for the uncoupled case, as it should). However, for 
< p, q < 1, F{\) as given by Eq. (|T2*|l yields in general 
a poor agreement with numerical simulations, as can be 
seen in Fig. ^ for different values of p. When A ~ 0, 
Eq. 1(12)1 predicts F ~ A/(l — 2p), which leads to obvi- 
ously nonphysical results for p > 1/2 (see leftmost part 
of Fig. In particular, F(X) is multi-valued, leading to 
lim A ^ + F^0. The mean field result therefore suggests 
a transition to an ordered state at A = which is simply 
forbidden by the automaton rules |27| . By generalizing 
Eq. i|ll|) , this failure to predict the absorbing state of the 
system can in fact be extended to regular lattices with co- 
ordination z, where the single-site approximation yields 



FIG. 2: Response curves for (a) n = 3 and (b) n — 10 au- 
tomata: simulations (symbols) and pair approximation [lines, 
according to Eqs. 1171181 1. From bottom to top, p = 0, 0.3, 
0.6 and 1, q = l-(l-p) 2 . The pair approximation eliminates 
the small- A anomalies of the single-site solution, yielding ex- 
cellent agreement with simulations for the extreme cases p — 
and p — 1. 



clearly not satisfactory for the calculation of the dynamic 
range, a refinement is needed. 



V. COUPLED NEURONS: PAIR 
APPROXIMATION 

The pair approximation consists in keeping the influ- 
ence of only one neighbor in the conditional probabili- 
ties 0, thus P t {j\\h, ■ ■ .,jm) ~ Pt(h\h)' In this case 
m-site probabilities are reduced to combinations of up to 
two-site probabilities. In particular, three- and four-site 
probabilities become pf| 



P{k,l,m) 
P(j,k,l,m) 



P(k,l)P(l,m) 

W) 

P(j,k)P(k,l)P(l,m) 



(13a) 
(13b) 



P(k)P(l) 

It is therefore possible to rewrite Eq. IJ^J in this approx 
imation: 



P t+1 (l) w AP t (0)+(l-A)P t (l,0) 



F 



A/(l — pz). Since this level of approximation is 



2p+{q-2p) 



#(1,0) 



Pt(0) 
(14) 
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Eq. (|14fl . on its turn, depends on Pt(l, 0), whose evolution 
can be exactly obtained (up to homogeneity and isotropy 
assumptions) : 

P t+1 (l,0) = \P t (n - 1,0) +p(l-\)P t (n- 1,0,1) 
+A(l-A)P t (0,0) 
+p(l- A)(l-2A)P t (l,0,0) 
-p 2 (l-A) 2 P t (l, 0,0,1). 

(15) 

With the help of the pair approximation in Eqs. <|13|) . 
Eq. I|15|) becomes 



P+i(l,0) w P t (n- 1,0) 



A + p(l - A) 



Pt(0) 



+(l-A)P t (0,0) 



A+p(l-2A) 



Pt(l,0) 
Pt(0) 



-p 2 (l-X) 



P t (l,0) 2 



P(0) 2 



(16) 



Since Pt(j, 0) depends on P t (j — 1, 0) and P t (j — l,n— 1), 
and Pt(0, 0) depends, among others, on P t (n — l,n — 1), 



all the equations for two-site probabilities are in princi- 
ple required for the dynamical description of the system. 
Together with the equations for single-site probabilities, 
they form a (n 2 + 3n)/2-dimensional map whose station- 
ary stable solution can be analytically studied. While 
the Appendix contains details of the derivation of those 
equations, we discuss the main results below. 



The main point to be noted is that the calculation of 
the stationary state presents additional difficulties when 
n > 4. In that case, the pair probabilities P(j, 0) with 
2 < j < n — 2 have the same stationary value, but differ 
from P(n — 1,0). In particular, for p = q = 1 one ob- 
tains P(j, 0) = [2 < j < n - 2, see Eq. (|A~10)l ]. which 
in turn leads to many other vanishing probabilities and 
gives the deterministic case a sparse stationary matrix 
[see Eqs. (TOll . (fOj) and (|A~gj) ]. Those terms do not 
exist for the n = 3 case, which makes its analysis consid- 
erably simpler. In either case, for n > 3 one obtains the 
reasonable result P(n — 1,0) « P(l, 0), the l.h.s. (r.h.s.) 
being associated to the end (beginning) of an excitable 
wave front [see Eq. l|A.12|) ]. Combining these results, a 
normalization condition and the linearity of Eq. <|16[) in 
P t (0,0), we obtain (see Appendix): 



P(0)-P(l,QH2 + (n-3) 



(l-p)P(0) + (p-g)P(l,0) 
P(0)-pP(l,0) 

P(l,0)P(0)[P(0)-pP(l,0)] 

AP(0) 2 + p(l - 2A)P(0)P(1, 0) - p 2 (l - A)P(1, 0) 2 ' 



(17) 



which is valid Vn > 3. Consider now the stationary state of Eqs. © and l|14f> . They can be combined in a quadratic 
equation for P(1,0), yielding 



(2p-g)P(l,0)wG±(P(0)) 



: pP{0) ± 



'P(O) {P(0)[(n - l)p 2 + 2p - q + A(n - l)(2p - p 2 - q)] + (q- 2p)} 



(n-l)(l-A) 
I 



(18) 



Since P(l, 0) must vanish Vp, q in the limit A 
the only acceptable solution. 



0, G_ is 



The solution of Eqs. l(T7|) and (fTH|l determines P(0) as a 
function of A. Instead of numerically solving them, we it- 
erate the (n 2 + 3n)/2-dimensional map involving the one- 
and two-site probabilities for each value of A until it con- 
verges to its stationary state. Despite the growing num- 
ber of equations with n, this method has the advantage 
of avoiding unstable fixed points 26] [Eqs. I|17|) and l|18l) 
can have more than one solution]. Once P(0) is known, 
the response P(l) = P(A) is obtained via Eq. ©. 



A. Deterministic spike propagation (p = 1) 

Ordinary differential equations (ODEs) are the stan- 
dard modeling tool in computational neuroscience. This 
is due to the fact that, despite the stochastic nature of 
the opening and closing of individual ionic channels, a 
neuron containing a large number of such channels can 
very often be extremely well described by a deterministic 
dynamics |l4[ (an approach which has been established 
since the seminal work of Hodgkin and Huxley (3). In 
the present context, it is therefore important to address 
the case p = 1. This limit is consistent with a vari- 
ety of scenarios in which, in addition to the dynamics 
of individual neurons, spike transmission is also well de- 
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FIG. 3: Linear- log plot of the response curve for n = 3 
automata with p = q = 1 (filled circles), p — 0.5 and 
q = 1 — (1 — p) 2 (open triangles), and p = q = (open 
circles). Lines correspond to the pair approximation. Hori- 
zontal lines are F = 0.1F max and F = 0.9F max , vertical lines 
are A = Ao.i and A = A0.9 and arrows illustrate the dynamic 
range 5\ [Eq. Q] for p — and p = 1. The dynamic range 
of a chain of neurons with deterministic spike propagation is 
about twice as large as that of its uncoupled counterpart. 
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FIG. 4: Dynamic ranges (triangles for S\, circles for S r ) as 
a function of the number of states of the GHCA, obtained 
from the stationary solution of the pair approximation. Open 
(filled) symbols correspond to the p = q = (p — q = 1) 
case. Inset: 5\ as a function of p for n — 10 for simulations 
(dashed line) and pair approximation (solid line). In spite 
of the underestimation of the response observed in Fig. 
the pair approximation is able to reproduce the behavior of 
the dynamic range as a function of the probability of spike 
transmission. 



scribed by deterministic behavior. Specifically regarding 
our present study, deterministic spike transmission due 
to electrical coupling has previously been employed in 
the literature to model axo-axonal interactions both via 
ephaptic interactions (e.g. in the olfactory nerve llll t 
and gap junctions (e.g. in the hippocampus [53, l29|). 



FIG. 5: Log-log plot of the response curve for p = q = 1. Pair 
approximation (solid lines) and simulations (symbols) follow a 
power law (a = 1/2) for weak stimuli, while finite size effects 
lead to a linear response F ~ LA (dotted lines) for A < A C (L). 



This is in contrast with, say, dendro-dendritic gap junc- 
tions or chemical synapses (in the latter case, synaptic 
transmission can sometimes be as low as 10% due to the 
inherent stochasticity in the process of neurotransmitter 
release jQ, H3)> where the p = 1 limit can hardly be 
expected to apply. As we shall see in the following, in 
addition to its biological relevance, the response func- 
tion for p — 1 also has a different characteristic exponent 
which will help us understand the limiting behavior for 

Figure |3 shows the excellent agreement between the 
pair approximation and the simulations when p = q = 1 . 
One observes that the response is particularly enhanced 
in the low stimulus range. This feature is best seen in 
the logarithmic scale of Fig. |3J in comparison with the 
uncoupled case p = 0, the effect of the electrical interac- 
tion is to increase the sensitivity of the response for more 
than a decade, leading to a dramatic rise of the dynamic 
range. 

For each value of n, we can thus obtain the station- 
ary response F(X) and the dynamic ranges 5\ and 8 r in 
the pair approximation. Even though the response curve 
changes considerably for varying n (since F is bounded 
by F max = 1/n, see Fig. the dynamic range levels 
off smoothly, as can be seen in Fig. 0] For increasing 
n, the dynamic range of the p = q = I case approaches 
twice the value for the uncoupled case. The fact that 
this result holds for both S r and 5\ can be understood on 
the basis of the low-stimulus amplification, which plays 
the central role in the phenomenon: in this regime A is 
approximately linear in r. Should one choose a different 
relationship A(r), S r would obviously have different val- 
ues, but the drastic enhancement in the response due to 
the electrical coupling would not be affected. 

In order to understand the low-stimulus amplification 
induced by the coupling, we have analyzed Eqs. (|17I18|) 
when A ~ 0. Inspection of Fig.|21and previous numerical 
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simulations 01 suggest that P(l) ~ CX a , with a < 
1. This ansatz can be inserted into Eqs. © and (flTt - 
I18JI for general p and g, yielding a = 1/2 and p = 1 
as solutions. Deterministic spike propagation therefore 
leads to a power law response 

F(X) A ~° V2A , (19) 

a result that holds Vn, <?, as should be expected. This 
power law suggests a Hill function with a = 1/2, which 
is an excellent approximation for F(X) in the whole A 
interval when n is large. This result explains the doubling 
of the dynamic range as compared to the uncoupled case 
and is reminiscent of reaction-diffusion processes modeled 
by lattice gases [101 03 an d partial differential 

equations |35|]. Since the Hill function can be regarded 
as a saturating Stevens law, it is interesting to note that 
the experimental values of the Stevens exponents for light 
and smell intensities are respectively a ~ 0.5 and a ~ 
0.6 0. 

Let us now consider a chain with finite L and a very 
small value of A such that a single external stimulus oc- 
curs in a given time interval. In this case, the determin- 
istic nature of the propagation would lead to L spikes in 
the chain, while a single spike would be observed if the 
neurons were uncoupled. One would thus have F ~ Lf, 

and since / ~ A [from Eq. ©] we obtain F ~ LA. 
This corresponds to a linear regime where excitable waves 
do not interact. If one increases A, waves will start anni- 
hilating each other, leading to the power law response 
of Eq. qi'jp. as can be clearly seen in Fig. [SJ For a 
given system size L, there is therefore a crossover value 
X C {L) ~ 2/L 2 from a linear to a nonlinear response. In 
an infinite chain, there is no linear response since for any 
nonzero stimulus rate two excitable waves will inevitably 
interact. 

To assess the finite size effects in the biological context 
of the model, we notice that the dynamic range will be 
affected only if X C (L) > Ao.i, that is, for L < 20n. For 
neurons with refractory periods of the order of tens of ms, 
neuronal assemblies with L > 10 3 ~ 4 should therefore be 
well approximated by the limit L — > oo, as can be seen in 
Fig. El It is important to emphasize, however, that even 
small chains dominated by finite size effects still possess 
dynamic ranges which are significantly larger than those 
of the uncoupled case. For Ao.i ^ X C (L), the dynamic 
range increases approximately logarithmically with the 
total number of connected neurons, a result which holds 
for regular lattices in any dimension 17]. 

B. Probabilistic spike propagation (p ^ 1) 

For p ^ 1 , communication between spiking and resting 
neurons may eventually fail. This provides us with the 
simplest test under which the robustness of the mech- 
anism for dynamic range enhancement can be checked. 
From the biological point of view, this regime could be 
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FIG. 6: Dynamic range as a function of the system size L for 
p — q — 1. Lines are just guides to the eye. 

useful for modeling networks of neurons connected by 
chemical synapses, for instance. 

We start the analysis of the p ^ 1 case by noticing in 
Figs.|21and|3]that the agreement between simulations and 
the pair approximation is better than the mean field re- 
sults (specially in the low-stimulus region), but certainly 
not so good as in the extreme cases p = and p = 1. This 
inevitably affects the estimation of the dynamic range 
via the stationary state of the pair approximation (see 
below), but nonetheless allows us to understand qualita- 
tively how the response changes as p varies. 

As pointed out in the preceding section, the dynamic 
range is enhanced for p — 1 primarily due to the low- 
stimulus amplification associated to the propagation of 
excitable waves. As opposed to the deterministic case, 
however, for p / 1 a single excitable wave traveling in 
an infinite chain initially at rest will eventually die out. 
We should therefore expect a qualitative change in the 
response function for A ~ 0. This is indeed confirmed by 
reinserting the ansatz P(l) ~ CX a in Eqs. JBJ and ifpTT - 
I18fl without the constraint a < 1 . In this case, the linear 
behavior suggested by the plots in Fig. is easily con- 
firmed: 

m >3°(l±l)x, ,20, 

which is again valid Vn, q. Therefore, the low-stimulus re- 
sponse for p < 1 is governed by a = 1 , which is confirmed 
by the simulations displayed in Figd Interestingly, such 
a change in exponent for p < 1 seems to be absent from 
reaction-diffusion models in lattice gases m m m ei 
as well as partial differential equations |35l ] . 

Thanks to the growing coefficient in Eq. 1)20(1 . for p < 1 
the proximity to the transition that occurs at p = 1 
produces a crossover in the response from a linear to 
a square root behavior, dismissing the suspicion that a 
larger exponent might severely deteriorate the enhance- 
ment of the dynamic range (see Fig. 0). In particular, 
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notice that, for a = l/2is the dominant exponent 

at F = 0.1F max , which is used to calculate the dynamic 
range (see horizontal arrow in Fig. UJ. This explains the 
smooth monotonic increase in S\ with p, as shown in the 
inset of Fig.^ even though the exponent changes discon- 
tinuously at p — 1. On the one hand, we observe that 
deterministic spike propagation (p = 1) is certainly not 
essential for the enhancement of the dynamic range, in 
the sense that any p > yields a better response than 
uncoupled neurons. On the other hand, it is interesting 
to point out that, as p is varied from to 1, the increase 
in dynamic range is particularly pronounced for p > 0.9. 
This is in agreement with the conjecture that the relia- 
bility of electrical coupling among spiking neurons could 
indeed play a significant role in early sensory processing. 

VI. CONCLUDING REMARKS 

We have calculated the collective response to a Pois- 
son stimulus of a chain of electrically coupled excitable 
neurons modeled by n-state Greenberg-Hastings cellular 
automata. The single-site mean field approximation has 
been shown to give poor results, failing to predict the 
absorbing state of the lattice in the absence of stimulus 
for p > 1/2. The pair approximation yields a response 
curve which agrees reasonably well with simulations in 
the whole stimulus range. It is interesting to remark that 
the agreement is particularly good when p = q — 1, a 
deterministic regime in which the GHCA lattice mimics 
a system of coupled nonlinear ODEs. This reinforces an 
interesting perspective in the context of computational 
neuroscience: the possibility of applying techniques from 
nonequilibrium statistical mechanics to the study of spa- 
tially extended nonlinear systems. 

The enhancement of the dynamic range in the presence 
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FIG. 7: Log-log plot of the response curve: pair approx- 
imation (solid lines) and simulations (symbols) with q = 
1 — (1 — p) 2 and n — 3. For p < 1, there is a crossover 
between a — 1 and a = 1/2. The horizontal arrow shows 

0. l-Pmax ■ 



of electrical coupling is due to low-stimulus amplification. 
For uncoupled neurons (p = 0) the response is governed 
by the Hill exponent a = 1, leading to a dynamic range 
of ~ 19dB. For coupled neurons this value can be dou- 
bled in the limit p = q — > 1, when the Hill exponent 
becomes a — 1/2. This value is close to Stevens ex- 
ponents observed in psychophysical experiments of smell 
and light intensities. For < p < 1, the exponent re- 
mains a = 1, but the dynamic range increases smoothly, 
which can be understood on the basis of the crossover 
behavior observed in the response function for p < 1 • 

In the context of experiments at the cellular level, the 
enhancement of the dynamic range associated with an 
increase in sensitivity is also observed in both the olfac- 
tory Q and visual Q systems. While the dynamic range 
of OSNs (the neurons which perform the initial trans- 
duction) is about ~ lOdB the glomeruli (the next 
processing layer) have dynamic ranges at least twice as 
large It remains to be investigated experimentally 
whether this enhancement is indeed due to ephaptic in- 
teractions among the unmyelinated OSN axons in the 
olfactory nerve. 

Stronger experimental support for our conjecture on 
the role of electrical interactions is available for the mam- 
malian retina. Deans et al. Q have measured the firing 
rates of on-center ganglion cells for varying light inten- 
sity (measured in isomerized molecules of rhodopsin per 
rod per second, or Rh*/rod/s). The response curves have 
been obtained for both wild type (WT) mice as well as 
mice in which the expression of the protein connexin36 
(responsible for the gap junction intercellular channels) 
has been genetically knocked out (Cx36-KO). The differ- 




r (Rh*/rod/s) 

FIG. 8: Experimental response curves (normalized firing rate 
vs. light intensity) of retinal on-center ganglion cells in linear- 
log (main plot) and log-log (inset) scales (data extracted from 
Fig. 6 of Ref. 0). Filled (open) circles are for WT (Cx36- 
KO) mice, solid (dashed) lines show the results of the pair 
approximation, thus L — > oo, with p = q = 1 (p = q = 0). 
Upper curves are for n = 10, lower curves are for n = 3. 
The dot-dashed line corresponds to simulations with n = 10, 
p = q — 1 and L — 20. 
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ence in the response curves can be seen in Fig. [S] They 
present the same qualitative behavior of the curves shown 
in Fig. |31 exhibiting an increase in dynamic range in the 
presence of electrical coupling: 14dB for Cx36-KO and 
23dB for WT, values which are of the same order as those 
of Fig. El In particular, the exponent of the "coupled" 
(WT) case is a ~ 0.58 (see inset), which is slightly larger 
than what is obtained in the pair approximation. 

The quantitative agreement between the analytical and 
experimental curves is limited. On the one hand, the 
theoretical n — 3 curve can provide a good fit of the 
Cx36-KO data for p = q = 0, while the coupled case 
p = q = 1 does not adjust well to the WT data. For 
n = 10 and p = q = 1, on the other hand, the WT data 
are well matched by simulations with a finite L = 20 
system (staying below the L — > oo pair approximation), 
but for p = q = the same n = 10 automata are unable 
to give a good fit of the Cx36-KO data. The difficulties 
of a quantitative match are not surprising: the retina 
is organized in layers which have, to first order, a two- 
dimensional structure, signal processing from the pho- 
toreceptors to the ganglion cells involves a complex in- 
termediate neuronal circuit (with bipolar, horizontal and 
amacrine cells [3(> ;. and individual neurons themselves 
can have subtle dynamical properties (such as adapta- 
tion, for instance). All these properties are clearly ab- 
sent from our simple one-dimensional CA model. Yet it 
correctly predicts the reduction in the dynamic range of 
a neuronal system which loses electrical coupling among 
its cells. 

In order to have a quantitative agreement between ex- 
perimental and theoretical curves, additional modeling 
efforts are needed which incorporate specific details of 
the system under consideration. However, the response 
of simple models of excitable media remains an impor- 
tant subject to be studied, precisely because they have 
the potential to reveal simple mechanisms and scaling re- 
lations [35[ whose robustness can thereafter be subjected 
to further testing in experiments and more detailed mod- 
els. In this context, the simple Greenberg-Hastings CA 
strikes an interesting balance, on the one hand capturing 
essential features of collective neuronal dynamics, while 
on the other hand lending itself to analytical techniques 
borrowed from nonequilibrium statistical mechanics. 



APPENDIX: THE EQUATIONS FOR TWO-SITE 
PROBABILITIES 



1. Dynamics 

In all derivations below, homogeneity and isotropy are 
assumed. The sign denotes that the equality holds in 
the pair approximation [Eqs. i|13|) ]. We start by writing 
down the equation for Pt(0, 0), which holds Vn > 3: 



P t +i(0,0) - P t (n-l,n-l) 

+2(1 - A) [P t {n - 1, 0) - pP t (l, 0, n - 1)] 
+ (1- A) 2 [P t (0,0)-2pP f (l,0,0) 
+p 2 P(l, 0,0,1)] 
« P t (n -l,n- 1) 

P*(1,0)" 



+2(l-A)P t (n-l,0) 



l-p 



+ (l-A) 2 P t (0,0) 
2 P(1,0) 2 ' 



1 - 2p 



P(0) 
P(1,0) 
P(0) 



-P 



P(0) 2 



(A.l) 



The dynamics for two-site probabilities in the refractory 
period obey a simple recursive rule due to the determin- 
istic evolution of the automata: 



P t+1 (j,k)=Pt(j-l,k-l) , 2<j,k<n-l. (A.2) 

On the one hand, diagonal terms Pt(J,j) with j > 2 
recursively depend on P t (l,l), whose dynamics can be 
written as follows: 



P t+1 (l,l) = A 2 p(0,0) + 2pA(l-A)P t (l,0,0) 
V(l-A) 2 P t (l, 0,0,1) 



Pt (0,0) 



A 2 + 2p\(l - A) 



P(1,0) 
P(0) 



+ p 2 (l-A) 



2 P(1,0) 2 



P(0) 2 



(A.3) 



Off-diagonal terms, on the other hand, ultimately depend 
on Pt(j, 1). For j = 2, the equation is simply 
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P+i(2,l) - (X+p-p\)P t (l,0) 

+ (1-A)( (? - P )P(1,0,1) 



P(l,0)<j A + (l-A) 
,P(1,0) 



+{q-p) : 



P(o) 



(A.4) 



while for j > 3 one has 
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p t+1 0'i) = \PtU- 1,0) +p(i-\)p t (j- 1,0,1 



PO'-i,o) 



A + p(l - A) : 



P(0) 



(A.5) 



Finally, one needs equations for P t (j, 0) , j > 2 [recall 
Eq. JXBJ for P t (l,0)]. Like in Eq. (jXtj) . the case j = 2 
must be considered separately: 

P t+1 (2,0) = P t (l,n-l) + (l-A)(l-p)P t (l,0) 
+(l-A)(p-?)P t (l,0,l) 



P t (l,n-l) + (l-A)P t (l,0) 
Pt(l,0)' 



+(p-<z) : 



P(o) 



(i-p) 

(A.6) 



For j > 3, on the other hand, one immediately obtains 



P t (i-l,n-l) 

+(1-A) [P t (3- 1,0) -pP t (j- 1,0,1)] 
P t (i-l,n-l) 

P(i,o) 



+(l-A)P t (j-l,0) 



1-p- 



P(0) 



(A.7) 



which completes the set of all pair equations. Upon it- 
eration of Eqs. H2l4ll4ll6IA.llA.7jl . normalization condi- 
tions properly imposed in the initial conditions are nat- 
urally preserved. To determine the response function 
P(l) = P(A), we wait until the (n 2 + 3n)/2-dimensional 
map reaches a stationary state for each value of A. We 
describe below how the analysis of the stationary state 
can be reduced to just two equations [Eqs. (|17I18|) ]. 



Notice that we have a nonhomogeneous set of n — 4 linear 
equations for Xj = P(j,0): Xj ax n -j + (1 — a)xj-i, 
where a = A + p(l - A)P(1,0)/P(0) and x 2 = P(2,0) 
accounts for the nonhomogeneity in the equations for 
X3 and x n -2- The solution of these equations is sim- 
ply x n -2 ~ x n -i ss . . . w X3 w X2, as can be checked by 
inspection. The combination of Eqs. (|A.6jl and l|A.5jl in 
the stationary state, on the other hand, leads to 



P(j,0) 



J[P(1,0),P(0)] 

p(i,o) f* 1 -^ 111 • 'A» /> < L,, » 



P(0)-pP(l,0) 



2 < j < n - 2 



(A.10) 



One therefore obtains 



P(n- 2,0) w P(n-3,0)w...wP(2,0) 

« J[P(1,0),P(0)] . (A.ll) 



Finally, notice that P(n —1,0) can be obtained by com- 
bination of Eqs. (|A.7jl . (|A.11|> and (|A.4jl : 



P(n-l,0)«P(l,0) 



(A.12) 



2. Stationary state 

We start by handling the case n > 4. In the stationary 
state, the first term on the r.h.s. of Eq. <|A. T|) becomes, 
via recursive iterations of Eq. (|A.2ll , 



(A. 



P(j -l,n-l) = P(l, 1 + n - j), Vj > 3 



The above result can on its turn be further developed by 
means of Eq. I|A.5|> as long as 1 + n — j > 3, rendering 
the stationary state of Eq. (| A. 7|> : 



P(j,0) « P(n-i,0) 



A + p(l - A) 



+(l-A)P(i-l,0) 
3 < j < n - 2 . 



1-p 



P(1,0) 
P(0) 

P(1,0) 



P(0) 



(A.9) 



which completes the proof for n > 4. For n = 4, it suffices 
to invoke Eqs. HA.6|) and (|A.5jl to show that P(2, 0) w 
J[P(1,0),P(0)]. With this result, Eq. ffKJ^i holds for 
n > 4. Finally, for n = 3, Eqs. HA.6|) and (|A.4|I together 
also lead to Eq. (|A.12(I . 

Invoking the normalization condition Pt(0) = 
^™ =0 Pt(j, 0), one can deduce that, on the one hand, 



P(0, 0) = P(0) - 2P(1, 0) - (n - 3) J [P(l, 0), P(0)] . 

(A.13) 

On the other hand, in the stationary state Eq. (|16|l de- 
pends linearly on P(0,0), so it can be inverted, yielding 
[after substitution of Eq. i|A.12|) ] P(0,0) as a function of 
P(1,0) and P(0). Equaling this function to Eq. (|A~T^jl . 
P(0, 0) is eliminated and one obtains Eq. (|T7jl . 
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